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ABSTRACT 

This paper focuses on the statistical properties of chaotic orbit ensembles evolved in 
triaxial generalisations of the Dehnen potential which have been proposed recently 
to model realistic ellipticals that have a strong density cusp and manifest significant 
deviations from axisymmetry. Allowance is made for a possible supermassive black 
hole, as well as low amplitude friction, noise, and periodic driving which can mimic 
irregularities associated with discreteness effects and/or an external environment. The 
chaos exhibited by these potentials is quantified by determining (1) how the relative 
number of chaotic orbits depends on the steepness of the cusp, as probed by 7, the 
power law exponent with which density diverges, and Mbh, the black hole mass; 
(2) how the size of the largest Lyapunov exponent varies with 7 and Mbh', and (3) 
the extent to which Arnold webs significantly impede phase space transport, both 
with and without perturbations. The most important conclusions dynamically are (1) 
that, in the absence of irregularities, chaotic orbits tend to be extremely 'sticky,' so 
that different pieces of the same chaotic orbit can behave very differently for times 
~ lOOOOt/) or more, but (2) that even very low amplitude perturbations can prove 
efficient in erasing many - albeit not all - of these differences. The implications of these 
facts are discussed both for the structure and evolution of real galaxies and for the 
possibility of constructing approximate near-equilibrium models using Schwarzschild's 
method. For example, when trying to use Schwarzschild's method to construct model 
galaxies containing significant numbers of chaotic orbits, it seems advantageous to 
build libraries with chaotic orbits evolved in the presence of low amplitude friction and 
noise, since such noisy orbits are more likely to represent reasonable approximations to 
time-independent building blocks. Much of the observed qualitative behaviour can be 
reproduced with a toy potential given as the sum of an anisotropic harmonic oscillator 
and a spherical Plummer potential, which suggests that the results may be generic. 
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1 MOTIVATION 

The work described here exploits recently developed ideas 
from chaos and nonlinear dynamics to better understand 
the dynamics of some seemingly realistic galactic potentials. 
These potentials reflect the fact that many/most early- type 
galaxies have a pronounced central density cusp (cf. Lauer et 
al. 1995), possibly associated with the presence of a super- 
massive black hole (cf. Kormendy & Richstone 1995); and 
that, at least for galaxies with comparatively shallow cusps, 
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there is often evidence for moderate deviations from axisym- 
metry (cf. Kormendy & Bender 1996). This work also em- 
braces the fact that real galaxies are continually subjected to 
various irregularities which the theorist might like to ignore, 
including 'high frequency' discreteness effects reflecting the 
existence of internal substructures and 'lower frequency' ef- 
fects reflecting, e.g., systematic pulsations or the effects of 
nearby objects. 

Recent interest in chaos in galactic dynamics has been 
driven primarily by data, both ground-based and from the 
Hubble Space Telescope (HST), which reveal that the den- 
sity of stars in early-type galaxies typically rises towards the 
center in a power-law cusp (Lauer et al. 1995, Byun et al. 
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1996, Gebhardt et al. 1996, Kormendy et al. 1996, Moller, 
Stiavelli, & Zeilinger 1995). For example, an analysis of more 
than 65 elliptical and SO galaxies has established that, at a 
resolution of < 0.1 arc-seconds, the surface brightness pro- 
file I(R) is best approximated by a power law profile oc i?~ 7 
with 7 ranging from near zero to unity. Most previous dy- 
namical studies of galaxies (e.g., the King models) assumed 
a constant density core, with a concomitant analytic central 
surface brightness, I(R) oc 1 - AR 2 H . The HST obser- 
vations require completely new dynamical models to predict 
kinematic properties of the central regions, and to ascertain 
whether supermassive black holes are actually present. 

There is also evidence that many galaxies may be more 
irregularly shaped than the nearly axisymmetric objects as- 
sumed as late as the 1970's. For example, twisted isophotes 
are interpreted as evidence for deviations from axisymmetry, 
and the existence of nontrivial residuals in a fit of the sur- 
face brightness distribution to a cos 29 law suggests further 
that many systems do not even exhibit the symmetries of a 
triaxial ellipsoid (cf. Bender et al. 1989). Indeed, these ob- 
served irregularities are so pronounced that they have been 
proposed as the basis of a new classification scheme for el- 
lipticals (Kormendy and Bender 1996). 

The crucial point, then, as stressed by Merritt and col- 
laborators (cf. Merritt 1996, Merritt & Fridman 1996), is 
that the combination of cusps and triaxiality seems to make 
chaos nearly unavoidable. In a non-cuspy triaxial galaxy, the 
central regions are dominated by regular box orbits with the 
topology of a three-dimensional Lissajous figure. Inserting a 
cusp or a supermassive black hole can destabilize these or- 
bits. One thus anticipates that many of the orbits passing 
close to the center of the galaxy must be chaotic, and that 
this feature could play an important role in the structure 
and evolution of these regions of the galaxy. Far from being 
something exotic and improbable, the vast majority of ellip- 
tical galaxies may contain large fractions of chaotic orbits. 

These observations run counter to the historical trend 
in galactic dynamics, the foundations of which go back over 
fifty years to a time when most astronomers had a physi- 
cal worldview that was dominated by integrable and near- 
integrable systems. In recent years, much attention has fo- 
cused on constructing self-consistent galactic models, ideal- 
ized as time-independent solutions to the collisionless Boltz- 
mann equation. In particular, given various simplifying as- 
sumptions about equilibrium shapes, one can construct ex- 
act analytic solutions such as the integrable Stackel (cf. de 
Zeeuw 1985) models. However, the new HST observations 
imply that the very idea of integrable self-consistent dynam- 
ical models must be rethought (cf. Merritt 1996). Indeed, as 
Gebhardt et al. (1996) put it, 'it seems very unlikely that ex- 
perience gained from the analysis of orbits in static Stackel 
potentials or of triaxial objects with analytic cores has much 
connection to the central regions of real galaxies.' 

Once it be admitted that stellar orbits in galaxies can 
have a large chaotic component, a host of fundamental 
questions arise that to date have had no fully satisfactory 
answers: How should one construct self-consistent equilib- 
ria with chaotic orbits? Are the fundamental time scales 
such as the Chandrasekhar relaxation time tR (cf. Chan- 
drasekhar 1943a) changed by chaos? On what time scale are 
the trapped chaotic orbits used (cf. Athanassoula et al. 1983, 
Wozniak 1993) to explain the shapes of certain galaxies un- 



stable? What is the effect of a large central point mass? 
How accurately can the invariant density associated with 
the chaotic part of the phase space be approximated? Will 
elliptical galaxies with cusps reach triaxial steady states or 
bypass them in favour of axisymmetric ones? These, and 
many other, questions are variations on the central theme: 
what are the dynamical consequences of chaos in galaxies? 

The objective here is to explore these issues for the tri- 
axial generalisations of the Dehnen (1993) potentials, which 
have been considered extensively by Merritt and collabora- 
tors (e.g., Merritt & Fridman 1996). These correspond to 
potentials generated self-consistently from the triaxial mass 
density 
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for c/a = 1/2 and (a 2 - b 2 )/(a 2 - c 2 ) = 1/2. Four cases 
will be discussed in detail, namely 7 = 0, 7 = 0.3, 7 = 1, 
and 7 = 2, ranging from no cusp to a very steep cusp, re- 
spectively. The analysis involves extracting the statistical 
properties of chaotic orbit ensembles evolved in this fixed po- 
tential both with and without low amplitude perturbations, 
including an analysis of what Merritt & Valluri (1996) have 
termed 'chaotic mixing' (cf. Kandrup & Mahon 1994a, Ma- 
hon et al. 1995, Kandrup 1998b). The aim is to understand 
both (i) the extent to which topological bottlenecks like an 
Arnold (1964) web can impede phase space transport in the 
unperturbed phase space and (ii) how even low amplitude 
irregularities can help orbits to traverse these bottlenecks. 
Assessing these topological effects is crucial for understand- 
ing the dynamics of potentials that admit a coexistence of 
both regular and chaotic orbits; and, as such, an important 
first step in determining whether it be reasonable to use 
them as viable candidates for collisionless equilibria. 

Basic questions to be addressed include the following: 

• To what extent does the efficiency of phase space transport 
depend on the steepness of the central cusp? In particular, 
does steepening the cusp, which appears to increase the over- 
all importance of chaos (cf. Merritt & Fridman 1996), also 
make phase space transport more efficient? 

• How is chaotic mixing impacted by low amplitude pertur- 
bations, modeled as friction and noise or periodic driving? 
In particular, how large do such perturbations have to be in 
order to have significant effects within a Hubble time tn 1 - 
Earlier work on generic complex potentials would suggest 
that even perturbations so weak as to be characterised by a 
relaxation time tR ~ 10 6 ir> or longer can be important dy- 
namically within lOOtu (Habib, Kandrup, & Mahon 1997). 

• How are bulk, statistical properties altered by the presence 
of a supermassive black hole? Earlier work suggests that a 
supermassive black hole can increase the overall abundance 
of chaotic orbits, and that it may render chaotic orbits more 
unstable, i.e., endow them with a larger Lyapunov exponent 
(cf. Udry & Pfenniger 1988, Hasan & Norman 1990), but it 
is not clear whether the presence of a black hole accelerates 
or impedes phase space transport. 

The answers to these questions will then be used to speculate 
about an even more fundamental issue, namely: does it seem 
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reasonable to expect that these, and similar, potentials can 
be used to construct self-consistent (near-)equilibria? 

Section 2 recalls some basic results from nonlinear dy- 
namics critical to a proper understanding of chaotic poten- 
tials with a complex phase space, neglecting the effects of a 
cusp and/or a supermassive black hole. Section 3 then de- 
scribes how, at least for the generalised Dehnen potentials, 
the insertion of a central cusp appears to alter the basic 
picture. Section 4 focuses on the stability of flows in these 
cuspy triaxial potentials towards low amplitude irregulari- 
ties. Section 5 then considers how the picture is complicated 
by the addition of a supermassive black hole. Section 6 con- 
cludes by summarising the principal results and speculating 
on their implications. 



2 PHASE SPACE TRANSPORT IN CHAOTIC 
HAMILTONIAN SYSTEMS 

The phase space associated with a Hamiltonian system ad- 
mitting only regular or only chaotic orbits tends to be com- 
paratively simple topologically. However, the coexistence 
of large measures of both regular and chaotic orbits leads 
generically to a complex phase space, the chaotic phase space 
regions being laced with a complicated pattern of cantori (cf. 
Percival 1979) or an Arnold (1964) web. The former arise 
for three-degree-of-freedom systems with two global isolat- 
ing integrals, e.g., axisymmetric configurations (as well as 
two-degree-of-freedom systems with only one isolating in- 
tegral); the latter for three-degree-of-freedom systems with 
only one isolating integral. 

It is often - but not always - true that, for fixed val- 
ues of the global isolating integrals, the chaotic phase space 
region is connected in the sense that a single chaotic orbit 
will eventually pass arbitrarily close to every point (strictly 
speaking [cf. Lichtenberg & Lieberman 1992], this neglects 
tiny chaotic regions nested inside KAM tori). However, in 
many cases one still finds that, over surprisingly long time 
scales, literally hundreds of dynamical times (to) or longer, 
the phase space is de facto divided into nearly disjoint re- 
gions by cantori or Arnold webs, which can serve as efficient 
bottlenecks. It follows that, over time scales of astronomi- 
cal interest, an ensemble of chaotic orbits, all with the same 
energy, can - even if there are no other isolating integrals - 
divide into two or more distinct populations, distinguished 
by (i) what parts of the chaotic phase space they occupy 
and/or (ii) how chaotic they are (cf. Mahon et al. 1995). 
This phenomenon can be probed and quantified using tools 
like short time Lyapunov exponents (cf. Kandrup & Mahon 
1994a), which were first introduced in a mathematically pre- 
cise setting by Grassberger et al. (1988) or, alternatively, by 
probing the extent to which the Fourier spectrum has broad 
band power (Kandrup, Eckstein, & Bradley 1997). 

Within a given phase space region unobstructed by ma- 
jor cantori or the Arnold web, chaotic mixing is usually quite 
efficient (Kandrup & Mahon 1994a, Merritt & Valluri 1996). 
In particular, for an initially localised ensemble of orbits 
one observes typically (i) an initial exponential divergence 
in phase space, followed by (ii) an exponential approach to- 
wards a near-uniform population of the phase space region, a 
so-called near-invariant distribution. Both these effects pro- 
ceed rapidly on a mixing time scale tm comparable to the 



natural time scale tL = 1/x associated with the largest Lya- 
punov exponent. However, cantori and Arnold webs can dra- 
matically increase the time required to approach a true in- 
variant distribution (Kandrup 1998b) . Separate phase space 
regions reach an 'equilibrium' comparatively quickly, but the 
time scale t eq for the different regions to communicate and 
reach a global equilibrium can be extremely long! 

The notion of an invariant distribution (cf. Lichtenberg 
& Lieberman 1992) , corresponding to a microcanonical equi- 
librium, i.e., a uniform population of the accessible phase 
space regions, plays a crucial role in any self-consistent equi- 
librium. When evolved into the future, a generic initial phase 
density /(0) will be transformed into a new f(t) 7^ /(0) 
and, as such, cannot serve as a time-independent build- 
ing block. However, a uniform phase density /i nv is invari- 
ant under an evolution governed by Hamilton's equations, 
i.e., /mv(i) = /inv(O). Such an /i nv thus constitutes a time- 
independent building block that can be used for constructing 
a self-consistent equilibrium. 

This is especially important for triaxial systems which 
typically admit only one global isolating integral (cf. Bin- 
ney & Tremaine 1987). Considering phase space distribu- 
tions that depend on one global integral and nothing else 
seems too restrictive to model centrally condensed triax- 
ial systems. (For example, any nonrotating equilibrium that 
depends only on energy must (cf. Perez and Aly 1996) be 
spherical.) This means that, if a triaxial potential admit- 
ting only one global integral is to be used to construct an 
equilibrium model, that equilibrium must involve 'nonstan- 
dard' time-independent building blocks, e.g., with different 
weights assigned to regular and chaotic orbits with the same 
energy (cf. Kandrup 1998a). The important point then is 
that these building blocks must be truly time- independent 
if the distribution is to correspond to a true equilibrium. 

One could envision near-invariant distributions / n i v cor- 
responding to nearly constant populations of chaotic phase 
space regions that are separated from the rest of the chaotic 
phase space by cantori or an Arnold web and, consequently, 
are nearly time-independent. However, as time elapses or- 
bits will diffuse through the surrounding obstructions and 
/niv will approach the true invariant distribution, /inv 

The possibility of distinguishing between 'sticky' and 
'not sticky' chaotic orbit segments was recognised nearly 
thirty years ago (cf. Contopoulos 1971). In particular, it 
was observed that a chaotic segment can be stuck near a 
regular island and behave as if it were regular for very long 
times, hundreds of dynamical times or longer, even though 
it will eventually become unstuck. Such near-regular chaotic 
orbits can play an important role in galactic modeling (cf. 
Athanassoula et al. 1983, Wozniak 1993, Kaufmann & Con- 
topoulos 1996). Conventional wisdom suggests that regular 
orbits must provide the skeleton for a self-consistent model 
(cf. Binney 1978). However, in certain critical phase space 
regions (e.g., near corotation) almost no regular orbits may 
exist, although sticky orbits are still abundant. Using near- 
regular sticky orbits seems a logical alternative. 

In principle this is completely reasonable, but there is an 
important potential concern: even very weak perturbations 
can dramatically accelerate phase space transport through 
cantori or along an Arnold web, allowing sticky orbits to 
become unsticky and vice versa (cf. Lieberman & Lichten- 
berg 1972, Lichtenberg & Wolf 1989, Habib, Kandrup, & 
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Mahon 1997). One obvious perturbation is discreteness ef- 
fects which (cf. Chandrasekhar 1943a), are usually modeled 
as dynamical friction and white noise (e.g., in the context 
of a Fokker-Planck description) . For regular orbits in a fixed 
potential, such friction and noise only have significant effects 
on a collisional relaxation time Ir which, for galaxies as a 
whole, is usually long compared with tn- However, such per- 
turbations can have significant effects on the statistical prop- 
erties of chaotic orbit ensembles on much shorter times by 
accelerating diffusion through cantori (Habib, Kandrup, & 
Mahon 1997) or along an Arnold web (Kandrup, Pogorelov, 
& Sideris 2000). 

Another class of perturbations reflects companion ob- 
jects/nearby galaxies and internal pulsations which, in at 
least some cases, can be modeled by a (near-) periodic driv- 
ing. Not surprisingly, such driving is most effective when the 
driving frequencies are comparable to the natural frequen- 
cies of the unperturbed orbits (cf. Kandrup, Abernathy, & 
Bradley 1995). However, driving can be important even if 
the driving frequencies are quite low compared with the nat- 
ural frequencies (cf. Lichtenberg & Lieberman 1992), the res- 
onant coupling in this case arising via subharmonics. Both 
these phenomena are known to be important in the physics 
of nonneutral plasmas (cf. Tennyson 1979, Rechester, Rosen- 
bluth, & White 1981, Habib & Ryne 1995). 

In a rich cluster, the external environment probably 
cannot be modeled simply as a superposition of a small num- 
ber of near-periodic forces. And similarly, there may be in- 
ternal irregularities that cannot be approximated reasonably 
as nearly periodic. However, one might still anticipate that 
these effects can be modeled as a 'random' (i.e., stochastic) 
process involving a superposition of a large number of dif- 
ferent frequencies (formally, any signal can be decomposed 
into a sum of periodic Fourier components). For a random 
combination of frequencies combined with random phases, 
this is equivalent mathematically to (in general) coloured 
noise (cf. van Kampen 1981, Honerkamp 1994), with a finite 
autocorrelation time and a band-limited power spectrum. 

It is significant that evidence for irregular shapes, as 
probed by isophotal twists or cos 30 corrections to isopho- 
tal ellipses, is especially common in high density environ- 
ments (cf. Zepf & Whitmore 1993, Mendes de Olivera & 
Hickson 1994). This suggests that such irregularities could 
be induced environmentally by collisions and other close en- 
counters; and, even more fundamentally, that these galaxies 
have been displaced from equilibrium by their surrounding 
environment. At the most extreme level, these observations 
could raise the question of whether it be reasonable to model 
galaxies as self-consistent equilibria. A more conservative re- 
sponse is to determine whether such self-consistent equilibria 
are stable towards low amplitude irregularities. 

The crucial point in all this is that even if, in the ab- 
sence of all perturbations, a galactic model behaves as a 
stable or near-stable equilibrium for times t > tH, it may be 
destabilised within t < tH by comparatively weak but 're- 
alistic' perturbations associated with internal irregularities, 
systematic pulsations, or the surrounding environment. 

The effects of these perturbations do not turn on 
abruptly: there is no obvious critical amplitude below which 
they are irrelevant. Instead, one finds that the efficacy of the 
perturbations typically scales roughly logarithmically with 
amplitude. Why do these perturbations have an effect? Pe- 



riodic driving is easily understood as inducing a resonant 
coupling between the driving frequencies and the natural fre- 
quencies of the unperturbed orbit. The driving is compara- 
tively ineffectual if these frequencies are extremely different, 
although one can get couplings via sub- and superharmon- 
ics. Noise-induced phase space transport can also be under- 
stood as a resonance phenomenon. Coloured noise with a 
band-limited power spectrum only has an appreciable effect 
if the noise has substantial power at frequencies comparable 
to the natural orbital frequencies. When the autocorrelation 
time t c is shorter than, or comparable to, to, so that there 
is appreciable power at frequencies ~ 1/to, coloured noise 
has almost the same effect as white noise. As t c increases to 
larger values, so that the power spectrum cuts off at lower 
frequencies, the efficacy of the noise decreases logarithmi- 
cally (Pogorelov & Kandrup 1999, Kandrup, Pogorelov, & 
Sideris 1999). 



3 UNPERTURBED DEHNEN POTENTIALS 

When considering orbits in a nonintegrable potential there 
are at least three distinct notions of 'how chaotic,' each of 
which is relevant to galactic dynamics. 

1) What fraction of the accessible phase space is chaotic 
and what fraction regular? Do there exist the requisite orbit 
families for the skeleton of a self-consistent model? Short 
of actually building a self-consistent distribution function 
there is in general no guarantee that any given potential can 
serve as an equilibrium - hence the importance of techniques 
like Schwarzschild's (1979) method (see also Schwarzschild 
1993). However, most dynamicists would probably agree 
that potentials for which the phase space contains only a 
very small measure of regular orbits are unlikely candidates 
for equilibria. 

2) How unstable are individual chaotic orbit segments, i.e., 
how large are the (short time) Lyapunov exponents? This 
question is important for at least two reasons: The size of 
these exponents regulates the rate at which nearby orbits 
diverge and, hence, the rate of chaotic mixing. Moreover, 
segments with especially small short time Lyapunov expo- 
nents may be nearly indistinguishable from regular orbits 
over times ~ tH, so that one might perhaps use them in lieu 
of regular orbits when modeling complex structures. 
These first two points are more or less obvious. A phase 
space can be 'very chaotic' in the sense that the Lyapunov 
exponents are very large but, nevertheless, 'not so chaotic' 
in the sense that the relative measure of chaotic orbits is 
comparatively small. Less trivial, perhaps, is the following: 

3) Are most of the chaotic phase space regions at fixed en- 
ergy connected in the sense that a single chaotic orbit can, 
and will, access the entire region? (Even neglecting tiny 
chaotic regions nested inside KAM tori, there is no guar- 
antee that this is so.) And, assuming that the answer is yes, 
to what extent can a single orbit pass unimpeded through- 
out that region without being obstructed by cantori or an 
Arnold web? Such questions related to phase space transport 
are important in galactic dynamics because they impact on 
the overall efficacy of chaotic mixing and the extent to which 
it makes sense to use 'sticky' chaotic orbits as near-regular 
building blocks. 
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3.1 What fraction of the orbits are chaotic? 

To obtain a reasonable estimate of the relative number of 
regular and chaotic orbits, one requires reasonable samplings 
of initial conditions. These were provided by generating or- 
bit libraries similar to, but more complete than, those com- 
puted by Merritt & Fridman (1996) for use in constructing 
self-consistent equilibria. This entailed approximating each 
model by 20 constant energy shells, corresponding to phase 
space hypersurfaces containing 1/21, 2/21, ... of the total 
mass. Attention focused primarily on three shells, namely 
the two lowest and the ninth lowest. The two lowest energy 
shells presumably feel the cusp most acutely; the ninth shell, 
corresponding to an intermediate energy, should be domi- 
nated less completely by the cusp. For each choice of 7 and 
energy E, orbits were generated for 1000 initial conditions, 
these corresponding closely to the classes of initial condi- 
tions originally considered by Schwarzschild (1979) (see also 
Schwarzschild 1993). 250 of the initial conditions had v = 
and, in the absence of chaos, would correspond presumably 
to box orbits. Another 250 initial conditions were chosen 
along each of the three principal planes, with the two com- 
ponents of velocity in the plane vanishing identically but the 
third component nonvanishing. These yielded orbits which, 
in the absence of chaos, might be expected to correspond to 
long, short, and intermediate axis tubes. 

Each orbit was integrated for a total time t = 200tu 
and an approximation to the largest Lyapunov exponent 
obtained by tracking simultaneously a small initial per- 
turbation that was renormalised periodically in the usual 
way (cf. Lichtenberg & Lieberman 1992). If the chaotic or- 
bits were not 'sticky,' one would expect such integrations 
to yield a relatively clean separation between regular and 
chaotic behaviour, integrations with other potentials indi- 
cating that, in this case, integrating for t = 100£d is usually 
adequate to distinguish between chaos and regularity (cf. 
Kandrup, Pogorelov, & Sideris 1999). However, for these 
triaxial Dehnen potentials such is not the case. In this case, 
segments of chaotic orbits, evaluated for times t = 200tn 
(or even much longer), exhibit invariably a broad range of 
short time Lyapunov exponents, extending from values of \ 
no larger than what is expected for regular orbits to sub- 
stantially larger values. For this reason, there can be some 
uncertainty in determining whether any given orbit is, or is 
not, chaotic. However, it is possible to determine the num- 
ber of "strongly chaotic" orbits with short time Lyapunov 
exponents larger than the near-zero values of \ computed 
for regular orbits. For the remainder of this subsection the 
appellation "chaotic" refers to such strongly chaotic orbits. 
The quoted values thus represent lower bounds on the actual 
number of chaotic orbits. 

For 7 = 0, there seem to be almost no (~ 12 out of 1000) 
chaotic orbits in the lowest energy shell. Approximately 10% 
of the orbits in the second shell are chaotic and the number 
increases to about 15% in the ninth shell. The near-absence 
of chaos in the lowest energy shell reflects the fact that, 
because there is no cusp, the very lowest energy orbits should 
be near-integrable boxes evolving in what is essentially an 
anisotropic harmonic potential. For 7 = 0.3, the relative 
fraction of chaotic orbits is quite similar except for the fact 
that, even in the lowest energy shell, about 5% of the orbits 
are chaotic. It is natural to attribute this small increase in 



the fraction of chaotic orbits to the fact that the smooth 
core has been replaced by a cusp, albeit a cusp sufficiently 
weak that the force acting on a star at r — > does not 
diverge. This interpretation is corroborated by the fact that, 
for 7 = 0.3, chaotic orbits in all three shells tend invariably 
to follow orbits that bring them close to the center. For 
example, the minimum distance from the origin for orbits 
in the lowest energy shell varies between r min = 0.003 and 
Train = 0.227, but all the seemingly chaotic orbits with \{t = 
200to) > 0.008 have r min < 0.07 and most have r m ; n < 
0.04. 

The same trend is also observed for the larger values of 
7, although the details are somewhat different. For 7 = 1, 
a significant fraction of the orbits is chaotic for all three 
energies but in this case the relative number of chaotic orbits 
decreases with increasing energy. Here approximately 25% 
of the lowest energy orbits are chaotic whereas the number 
decreases to ~ 15% for orbits in the ninth shell, roughly the 
same as for 7 = and 0.3. For 7 = 2, the relative fraction of 
chaotic orbits decreases from about 30% in the lowest energy 
shell to 25% in the ninth shell. This decrease presumably 
reflects the fact that, for 7 = 1 and 7 = 2, the central cusp 
plays a dominant role in generating chaos. (Recall that the 
force diverges at r — > for 7 > 1.) 

To say that much of the chaos is triggered by a 'close 
encounter' with the central cusp seems an oversimplifica- 
tion. For all three nonzero values of 7, there exist significant 
numbers of orbits with small r m i n that are unquestionably 
regular. Moreover, even for the cuspless model with 7 = 0, 
the chaotic orbits tend to have relatively small values of 
r m in- More reasonable is Merritt's (1996) argument that, at 
least for the lower energy shells, much of the chaos is trig- 
gered by an appropriate resonance overlap, which is stronger 
for steeper cusps. 

The relative numbers of orbits with different values of 
X can be gauged from FIGURES 13 and 14, which will be 
discussed more carefully in Section 5. 

3.2 How unstable are these chaotic orbits? 

To appreciate the significance of this chaotic behaviour, it 
is important also to assess the natural time scale associated 
with the exponential instability of these chaotic orbits, both 
in absolute units and expressed dimensionlessly in units of 
the dynamical time, to- The former provides information 
about the time scale on which the effects of chaos could be 
manifested physically, e.g., in the context of chaotic mixing. 
The latter ties the chaos more securely to the underlying 
dynamics. An estimate of the true Lyapunov exponent, de- 
fined in an asymptotic t — > 00 limit, was obtained for each 
value of 7 and E by computing \ f° r each of 8 chaotic ini- 
tial conditions with given 7 and E, integrated for a total 
time t — 102400to, and then constructing the average value 
(x), as well as dimensional and dimensionless time scales 
tL — l/(x) ana - Tl = 1/((x)*-d)- These integrations were 
performed using a variable time-step integrator which typi- 
cally conserved energy to at least one part in 10 5 or better. 

The principal conclusion of these computations is that, 
even though the dimensional Lyapunov time tL varies enor- 
mously for different choices of 7 and E, the dimensionless Tl 
does not. For the twelve different samples - four choices of 
7 and three choices of E - the computed values of Tl varied 
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Figure 1. (a) JV[x(At)], the distribution of short time Lyapunov 
exponents, for chaotic orbits in the lowest energy shell with 7 = 
for sampling time At = lOOto- The error bars reflect formal N 1 / 2 
uncertainties in the number of orbits in each bin. (b) 7V[x(At)] for 
the same orbits for At = 800t£>. (c) JV[x(Ai)] for chaotic orbits in 
the ninth energy shell with 7 = and At = 100t D . (d) 7V[x(At)] 
for the same orbits for At = 800fo. (c) - (h) The same as (a) - 
(d) for 7 = 0.3. 

by less than a factor of 2.5. In every case the Lyapunov time 
on which small perturbations grow exponentially is of order 
3 — Ito- This is consistent with work in other potentials, 
in both two and three dimensions, where, for systems ex- 
hibiting global stochasticity, the Lyapunov time is typically 
comparable to, but somewhat longer than, a characteristic 
crossing time (cf. Mahon et al. 1995, Kandrup, Pogorelov, & 
Sideris 1999). This suggests, however, that chaos should be 
much more important physically at low energies and/or for 
steep cusps, where the dynamical time to is comparatively 
short. Assuming, e.g., that the mass M — 5 x 1O 11 M and 
a — 5 kpc, one finds (cf. eq. [13] in Merritt & Fridman 1996) 
that, for 7 = 2, the dynamical time for the lowest energy 
shell is as short as 1.2 x 10 6 yr, whereas, for the ninth energy 
shell for 7 = 1, t D = 8.8 x 10 7 yr. 

3.3 Phase space transport 

Visual inspection of the long time trajectories of the chaotic 
orbits described above suggests that (most of) the chaotic 
portions of the phase space at any given energy are con- 
nected in the sense that a single orbit can and (presumably) 
eventually will pass arbitrarily close to every point in the re- 



gion. However, visual inspection of these trajectories also in- 
dicates that, for all four values of 7, but especially the larger 
values, chaotic orbits can be extremely 'sticky', even when 
viewed over intervals as long as lOOOOtc or more. Sometimes 
the orbit segment will closely resemble a near-regular box or 
tube; sometimes it will be wildly chaotic. These differences 
are manifested in the behaviour of short time Lyapunov ex- 
ponents, which can differ significantly for very long times. As 
in other potentials (cf. Kandrup, Eckstein, & Bradley 1997), 
one discovers that segments that look nearly regular tend 
to have relatively small short time x' s > whereas segments 
that are manifestly irregular tend to have larger values of x- 
The net result is that, for a broad range of sampling inter- 
vals At, the distribution of short time Lyapunov exponents, 
JV[x(At)], generated from the aforementioned long time in- 
tegrations, is distinctly bimodal. 

As is evident from Section 2, this fact is not surpris- 
ing. However, there is at least one important respect in 
which chaotic segments computed in these Dehnen poten- 
tials differ from, e.g., orbits in the generalised dihedral po- 
tential explored by Kandrup, Pogorelov, & Sideris (1999). In 
that potential, as well as many other generic potentials, one 
finds that, for intervals At > lOOOtu or so, the distinctions 
between 'sticky' and 'wildly chaotic' tend to be erased, so 
that different chaotic segments with the same energy have 
similar statistical properties. In particular, even if for small 
At the distribution N[x(At)] is bimodal, for At > 1000t D 
the distribution tends to be singly peaked. For these tri- 
axial Dehnen potentials this is no longer the case. Even 
At = 20000tD is not long enough for these distinctions to 
be erased. 

The fact that A r [x(At)] can be bimodal for compara- 
tively long times is illustrated in FIGURES 1 and 2 which, 
for a variety of different values of 7 and E, exhibits N[x(At)] 
for both At = lOOtc and At = 800t D . 

For significantly larger values of At, it becomes pro- 
hibitively expensive computationally to compute enough or- 
bit segments to generate a meaningful distribution of short 
time x' s - However, it is still possible to extract useful in- 
formation about the distribution by determining how the 
dispersion a x varies as a function of At. As discussed more 
carefully in Kandrup, Pogorelov, & Sideris (1999) (see FIG. 
3 in that paper, as well as Kandrup & Mahon 1994b), an ap- 
plication of the Central Limits Theorem (cf. Chandrasekhar 
1943b, van Kampen 1981) suggests that if, over the time 
scales of interest, chaotic orbit segments constitute a single 
population with the overall instability of the orbit at any 
two instants essentially uncorrelated, <r x should decrease as 
At -1 / 2 , whereas the existence of multiple populations im- 
plies a dispersion that decreases more slowly. The principal 
conclusion here is that, for chaotic orbits in these gener- 
alised Dehnen potentials, a x typically decreases very slowly, 
much more slowly than At~ 1/2 . For At < 2 5 t D = 32t D or 
so, a x typically decreases in a fashion roughly consistent 
with a At~ 1/2 dependence, but for 2 5 t D < t < 2 17 t D , i.e., 
32tu < t < 65536tu, a x decreases much slower. One thus 
anticipates that, consistent with the fact that, even for pe- 
riods 3> lOOOOto, different segments exhibit significant vari- 
ability in their visual appearance, there exist distinct pop- 
ulations that persist for t > 10000t D . FIGURES 3 and 4 
exhibit a x (At) for a variety of different values of 7 and E. 
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Figure 3. (a) The dispersion a x for chaotic orbits in the lowest 
energy shell with 7 = 0, computed as a function of sampling 
time At. (b) a x for chaotic orbits in the lowest energy shell with 
7 = 0.3. (c) The second energy shell with 7 = 0. (d) The second 
energy shell with 7 = 0.3. (e) The ninth energy shell with 7 = 0. 
(f) The ninth energy shell with 7 = 0.3. 



Figure 2. (a) N[x(At)], the distribution of short time Lyapunov 
exponents, for chaotic orbits in the lowest energy shell with 7 = 1 
for sampling time At = 100t£>. (b) N[x(At)] for the same orbits 
for At = 800*15- ( c ) N[x(At)] for chaotic in orbits in the ninth 
energy shell with 7 = 1 and At = 100t D . (d) N[x{At)] for the 
same orbits for At = 800fij. (c) - (h) The same as (a) - (d) for 
7 = 2. 



3.4 Chaotic mixing 

These results have obvious implications for chaotic mixing. 
Earlier work (cf. Kandrup & Mahon 1994a, Merritt & Valluri 
1998, Kandrup 1998b) would suggest (i) that, as probed by 
quantities like the phase space dispersions a x and a px , ini- 
tially localised ensembles of chaotic orbits will diverge expo- 
nentially, and (ii) that, as probed by various lower-order mo- 
ments and/or coarse-grained reduced distribution functions, 
they will exhibit a roughly exponential evolution towards 
a near-equilibrium (a near-invariant distribution), both ef- 
fects proceedings on a mixing time scale tu that correlates 
with the Lyapunov time scale tz, = 1/x associated with the 
largest Lyapunov exponent. However, to the extent that the 
typical values of short time Lyapunov exponents can differ 
substantially for different ensembles, one might anticipate 
that these ensembles could approach a near-equilibrium at 
very different rates; and, to the extent that orbits can be 
'stuck' in a given portion of the chaotic phase space for a very 
long time, one would not expect that the near-equilibrium 
should be independent of the choice of chaotic ensemble. 
Rather, one would anticipate a two- (or more) stage evolu- 
tion, whereby a rapid approach towards a uniform popula- 
tion of the easily accessible regions is followed by a slower 
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Figure 4. (a) The dispersion a x for chaotic orbits in the lowest 
energy shell with 7 = 1, computed as a function of sampling time 
At. (b) <r x for chaotic orbits in the lowest energy shell with 7 = 2. 
(c) The second energy shell with 7 = 1. (d) The second energy 
shell with 7 = 2. (e) The ninth energy shell with 7 = 1. (f) The 
ninth energy shell with 7 = 2. 
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Figure 5. x — y coordinates for a 1600-orbit ensemble evolved 
with 7=1 with the lowest energy recorded at six different times, 
(a) t = t D . (b) t = 4t D . (c) t = 8t D . (d) t = 16t D . (c) t = 40t D . 
(f) t = W0t D . 
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Figure 6. y — z coordinates for the same 1600-orbit ensemble 
recorded at six different times, (a) t = to- (b) t = 4trj. (c) t = 
8t D . (d) t = 16t D . (c) t = 40t D . (f) t = 10Ot D . 



evolution as orbits diffuse along the Arnold web to probe the 
entire accessible phase space. The approach towards a near- 
invariant distribution typically proceeds on a comparatively 
short time scale tm ~ 1/x- The final approach towards a 
true invariant distribution often proceeds on a time scale 

These ideas were tested for the triaxial Dehnen poten- 
tials by performing simulations which tracked the evolu- 
tion of ensembles of 1600 initial conditions localised within 
phase space cells of typical size 0.01 or less, and these ex- 
pectations were all confirmed. Numerical evolution of lo- 
calised ensembles of initial conditions indicates that quanti- 
ties like the dispersion a x in the coordinate x do indeed di- 
verge exponentially at a rate that correlates with an average 
short time Lyapunov exponent for the ensemble. Moreover, 
coarse-grained reduced distribution functions constructed 
from these orbits appear to converge exponentially towards 
a nearly constant /i nv . However, the rate at which the dis- 
persions grow and the rate of convergence towards /i nv can 
both depend on the choice of ensemble as well as the di- 
rection^) in phase space that are probed. For example, for 
the lowest energy shell with 7 = 1, there is often a propen- 
sity for ensembles starting with small z to diverge compar- 
atively slowly in the z-direction, so that a z approaches a 
near-constant value more slowly than a x and a y , and the 
coarse-grained distribution function converges towards /i nv 
more slowly in the z- and t> z -directions than in any other 
direction. This is illustrated in FIGURES 5 and 6, which 
exhibit, respectively, the x — y and y — z coordinates for the 
same 1600 orbit ensembles at times t = to, 4to, 8to, 16irj>, 
40ijj, and lOOto- Alternatively, the spatial dispersions a x 



and a z are exhibited in the top two panels of FIGURE 8, 
which is discussed more carefully in the following Section. 

That the near-invariant distributions generated from 
different chaotic ensembles with essentially the same energy 
need not be statistically identical is illustrated by the fact 
that the moments associated with different ensembles need 
not evolve towards the same near-constant values. An ex- 
ample thereof is exhibited in the top six panels of FIGURE 
7, which compares the time-dependent dispersions a x , a y , 
and cr z for two different ensembles of 6400 orbits, each sam- 
pling the lowest energy shell with 7 = 1. Both ensembles 
were generated from initial conditions with v = 0, i.e., or- 
bits which, in the absence of a cusp, might be expected to 
correspond to regular boxes. Each panel is indicative of a 
moment converging towards a near-constant value, but it is 
evident that these values are distinctly different for the two 
different ensembles. One might perhaps worry that these dif- 
ferences could reflect the fact that some of the orbits in these 
ensembles are actually regular. That this is not the case is 
clear from the last two panels in FIGURE 7, which exhibit 
distributions of short time Lyapunov exponents, N[x(t)], for 
each of the 6400-orbit ensembles, computed at t — 200to- 
In each broad range of values are assumed, but the 

distributions are unimodal and the lowest x is significantly 
displaced from \ = 0. 

4 THE EFFECTS OF INTERNAL AND 
EXTERNAL IRREGULARITIES 

The objective here is to determine how the statistical prop- 
erties of chaotic orbits evolved in these generalised Dehnen 
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Figure 7. (a) The dispersion a x (t) generated for one 6400 orbit 
chaotic ensemble sampling the lowest energy shell with 7 = 1. 
(b) tr x for another chaotic ensemble with the same E and 7. (c) 
and (d) u y for the same ensembles, (e) and (f) a z for the same 
ensembles, (e) A distribution of largest short time Lyapunov ex- 
ponents, N[x(t)], for t = 200to, generated for the first ensemble, 
(f) The same for the second ensemble. 



potentials change if the evolution equations are modified 
to allow for low amplitude irregularities which could mimic 
the effects of internal substructures and/or an external envi- 
ronment. This entailed repeating the simulations of chaotic 
mixing for unperturbed systems described in the preceding 
Section but now allowing for low amplitude perturbations. 
Three specific effects were considered: 

• Periodic driving, intended to mimic the effects of one or 
two companion objects and/or internal pulsations. This in- 
volved allowing for an evolution equation of the form 



dt 2 



av(r) 

dx a 



+ A a sm(u)at + <p a ), (a = x,y,z), (3) 



incorporating simple pulsations in three orthogonal direc- 
tions with different amplitudes, frequencies, and phases. 
• Friction and additive white noise, intended to mimic 
discreteness effects, i.e., gravitational Rutherford scatter- 
ing. This involved solving a Langevin equation (cf. Chan- 
drasekhar 1943b, van Kampen 1981) 



d 2 x a 



dV(r) 



T)V + F , 



(a = x,y,z), 



(4) 



dt 2 dx a 

with n a constant coefficient of dynamical friction and F 
a 'stochastic force.' Assuming in the usual fashion that F 
corresponds to homogeneous Gaussian noise, its statistical 



properties are characterised completely by its first two mo- 
ments, which take the form 

(F a (t)} = and 

(F a (ti)F b (t2))=S ab K(ti-t 2 ), (a,b = x,y,z). (5) 

The assumption that the noise be white implies that the 
autocorrelation function K is proportional to a Dirac delta, 
so that 



K(t) = 2t]@6 d (t). 



(6) 



The normalisation here ensures that the friction and noise 
are related by a Fluctuation-Dissipation Theorem at a tem- 
perature O. 

• Coloured noise, intended primarily to mimic the effects 
of a stochastic external environment. This entailed allowing 
for forces that are random but, nevertheless, of finite du- 
ration, i.e., characterised by a finite autocorrelation time. 
The specific example considered here corresponds to the so- 
called Ornstein-Uhlenbeck process (cf. van Kampen 1981), 
for which K decreases exponentially in time, i.e., 



K(r) = anO exp(— a|r|), 



(7) 



where the autocorrelation time t c = 1/a sets the time scale 
on which F changes appreciably. The normalisation here en- 
sures that the diffusion constant D that would enter into a 
Fokkcr-Planck description, 



j — < 



drK(j) = 2®n, 



(8) 



is independent of a. 

White noise was implemented computationally using an 
algorithm described in Griner, Strittmatter, & Honerkamp 
(1988). Coloured noise was implemented using an algorithm 
similar to that described in Pogorelov & Kandrup (1999). 

This modeling of perturbations is clearly simplistic but, 
nevertheless, not completely unreasonable. Earlier work on 
periodic driving as a source of phase space transport sug- 
gests strongly that the exact way in which the system is 
pulsed matters less than the amplitude and driving fre- 
quency. For example, Kandrup, Abernathy, and Bradley 
(1995) found that, for variants of an anisotropic Plummer 
potential, similar results obtained when pulsing the core ra- 
dius and the anisotropy parameter. 

Although the detailed forms of the friction and noise 
can be very important for phenomena like evaporation from 
a cluster (cf. Chandrasekhar 1943a) or barrier penetration 
problems (cf. Lindenberg & Seshadri 1981), they appear to 
be comparatively unimportant in determining the rate of 
phase space transport. Work on simpler two- and three- 
dimensional systems (Pogorelov & Kandrup 1999, Kandrup, 
Pogorelov, & Sideris 1999) shows that making n a simple 
but nontrivial function of v and/or turning off the friction 
term altogether has only minimal effects. Moreover, the de- 
tailed form of the coloured noise seems largely immaterial. 
What really matter are the amplitude and the autocorrela- 
tion time t c , which determines the frequencies at which the 
perturbation has substantial power. When t c < to, colour 
matters very little and the effects are similar to white noise. 
However, when t c 3> to, so that there is little power at fre- 
quencies ~ 1/to, the efficacy of the perturbation decreases 
logarithmically. 
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Figure 8. (a) The dispersion u x computed for an initially lo- 
calised ensemble of 1600 orbits in the lowest energy shell with 
7 = 1 integrated in the absence of any perturbations, (b) The 
corresponding a z . (c) a x for the same ensemble, now allowing for 
friction and white noise with temperature = — E and amplitude 
r) = 10~ 7 . (d) a z for the same friction and noise, (e) a x for the 
same ensemble, with = — E and r) = 10~ 6 . (f) u z for the same 
friction and noise, (g) u x for the same ensemble, with = — E 
and 77 = 10 -5 . (h) a z for the same friction and noise. 



4.1 Friction and white noise 

Consider first discreteness effects, modeled as friction and 
white noise. Here the most obvious conclusion is that even 
weak perturbations can significantly accelerate phase space 
transport by allowing orbits that are comparatively 'sticky' 
in one or more phase space directions to become 'unstuck.' 
This is illustrated in FIGURE 8, which exhibits the configu- 
ration space dispersions a x and a z for the ensemble of initial 
conditions used to generate FIGURES 5 and 6, contrast- 
ing unperturbed motion with orbits perturbed by friction 
and additive white noise with constant O = — E and three 
different values of 77, namely 77 = 0.4 x 1CT 7 , 0.4 x 10~ 6 , 
and 0.4 x 10~ 5 . Given that, for this choice of 7 and energy, 
to ~ 2.444, these values correspond respectively to relax- 
ation times ( H = l/i)W W 7 t D , l(ft D and 10 5 t D . 

It is apparent that friction and noise with 77 = 0.4 x 
10~ 7 has a relatively small effect, that 77 = 0.4 x 10 -6 has 
a substantially larger effect, and that perturbations with 
77 = 0.4 x 10~ 5 suffice to make a z evolve as if the ensemble 
were hardly sticky at all. This behaviour is also manifested 
by the visual appearance of the orbit ensemble, as illustrated 



-0.24 



n 0.00 



n 0.00 



-0.24 




0.24 



0.00 



-0.24 



(b) 



- n 0.00 




0.24 



n 0.00 



-0.24 




-0.24 0.00 0.24 

y 

Figure 9. y — z coordinates for the 1600-orbit ensemble of FIGS. 

5 and 6, now evolved with friction and white noise with = —E 
and X] = 10~ 6 . (a) t = t D . (b) t = 4t D . (c) t = 8t D . (d) t = 16t D . 
(c) t = 40t D . (f) t = lOOtD. 

in FIGURE 9, which, for 77 = 10~ 6 , exhibits the y and z 
coordinates of the particles at the same times as for FIGURE 
6. 

Even in the absence of stickiness, discreteness effects 
can be important by serving to 'fuzz out' high frequency 
structure in quantities like the phase space dispersions and 
other lower order moments, as well as various reduced distri- 
bution functions, thus rendering more efficient the approach 
towards a near-equilibrium. In particular, such noise can 
serve to damp the oscillations associated oftentimes with 
the evolution towards equilibrium, as noted, e.g., by Merritt 

6 Valluri (1996) and Kandrup (1998b). This effect is some- 
what apparent in FIGURES 8, which exhibit a x and o z for 
a period t — lOOto, but is even more striking in FIGURES 
10, which exhibit a y for the shorter interval < t < 50ti). 

How large must 77 be to have a significant effect? In 
general, perturbations corresponding to tR ~ 10 7 to tend to 
have (at least) a noticeable effect within ~ 20 — 30to, and 
perturbations with tn ~ 10 to can have an appreciable ef- 
fect on the overall approach towards a near-equilibrium on 
a comparable time scale. Perturbations with tR as small as 
can have drastic effects within 10to or so. 

Finally, consistent with Habib, Kandrup, & Mahon 
(1997), it was found that, for chaotic ensembles, the root 
mean squared change in energy scales as 

8E rms m (|£|e77t) 1/2 , (9) 

with E the original energy, which implies that noise is im- 
portant already for 8E rms ~ 10~ \E\. This reinforces the 
interpretation (cf. Pogorelov & Kandrup 1999, Kandrup, 
Pogorelov, and Sideris 1999) that this phase space transport 
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Figure 10. (a) The dispersion c y computed for an initially lo- 
calised ensemble of 1600 orbits in the lowest energy shell with 
7 = 1 integrated in the absence of any perturbations, (b) <T y 
for the same ensemble, now allowing for friction and white noise 
with temperature © = —E and amplitude r\ = 10~ 7 . (c) a y for 
the same ensemble, with = — E and r\ = 10~ 6 . (d) a y for the 
same ensemble, with = — E and n = 10~ 5 



is not driven by changes in energy, which only become large 
on a relaxation time tn, but instead reflects diffusion along 
Arnold webs on (nearly) constant energy hypersurfaces. 

The basic conclusion derived from these investigations 
is that discreteness effects can be important within a period 
as short as 50— lOOi o provided only that the relaxation time 
is no larger than 10 6 — 10 7 tc or so. This suggests strongly 
that, especially near the centers of cuspy galaxies, discrete- 
ness effects should be important on time scales considerably 
shorter than tn, the age of the Universe. 



4.2 Periodic driving 

Overall, periodic driving has the same qualitative effects on 
orbit ensembles as does white noise. Periodic driving can 
enable "sticky" orbits to become unstuck, thus facilitating 
phase space transport; and, even in the absence of topolog- 
ical obstructions, driving can help fuzz out high frequency, 
short wavelength structures. Not surprisingly, for fixed am- 
plitude A periodic driving has the largest effect for driving 
frequencies uji with ~ to, so that there is a substantial 
resonant coupling between the driving frequencies and the 
natural frequencies of the unperturbed orbits. 

One important difference between noise and periodic 
driving is that, whereas noise represents an incoherent pro- 
cess, involving a random superposition of frequencies, pe- 
riodic driving with fixed frequencies is a coherent process. 
That periodic driving is coherent would suggest that its cou- 
pling to orbits should scale linearly in amplitude, i.e., oc A, 
a fact that has been confirmed numerically (cf. Kandrup, 
Abernathy, & Bradley 1995). Because noise is incoherent, its 
coupling scales instead as n 1 ^ 2 (cf. Habib, Kandrup, & Ma- 
hon 1997). For 1/uji ~ to it is thus natural to expect that, 
if noise requires an amplitude r\ ~ 1CP P to induce a sub- 
stantial effect within a given time interval, periodic driving 
will require an amplitude A ~ 10~ p / 2 to have a comparable 
effect. 



Overall, this scaling tends to hold at least roughly true. 
However, one discovers (i) that the minimum amplitude re- 
quired to achieve a given effect can exhibit a comparatively 
sensitive dependence on the choice of frequencies, and (ii) 
that the "typical" effect of periodic driving is slightly weaker 
than is predicted by this scaling. These facts are both easily 
understood. The fact that noise involves a superposition of 
many different frequencies implies that the details tend to 
"wash out," and that one is more likely to have power at 
one of the "more effective" frequencies. 



4.3 Friction and coloured noise 

As for other two- and three-dimensional potentials, coloured 
noise sampling the Ornstein-Uhlenbeck process has virtually 
the same effect as does white noise with the same amplitude, 
provided only that the autocorrelation time t c is short com- 
pared with the dynamical time to- However, when t c be- 
comes comparable to to the effects of the noise begin to de- 
crease significantly; and, for t c ^> to the effects of the noise 
are substantially reduced. 

Examples of this behaviour, which corroborate the in- 
terpretation of noise-induced phase space transport as a res- 
onance phenomenon, are presented in FIGURE 11, which 
exhibits <j z generated from the same ensemble of initial 
conditions used to generate FIGURE 8, now evolved with 
coloured noise. The first three panels represent friction and 
noise with O = E and r\ — 0.4 x 10~ 5 for three choices of 
autocorrelation time, namely t c ~ QAto, \.2to, and 4.0to- 
The fourth and fifth panels were generated for O = —E and 
t c w 4.0tc but larger values of rj = 0.4 x 10~ 4 and 0.4 x 10~ 3 . 
The final panel is generated for 6 = — E, r) = 1.2 x 10~ 3 , 
and t c x Ylto- It is clear that for these values, which are not 
inappropriate for modeling the effects of an external environ- 
ment on a galaxy embedded in a rich cluster, perturbations 
can have a dramatic effect. 



4.4 The approach towards an invariant 
distribution 

Even though low amplitude perturbations can accelerate the 
approach towards a near-invariant distribution, they need 
not suffice to make different ensembles evolve towards the 
same distribution within a time t ~ lOOtc or so. For ex- 
ample, even though friction and noise corresponding to a 
relaxation time as long as ta ~ lQ 6 to or more can have an 
appreciable effect on the efficacy with which the ensemble 
approaches a near-invariant distribution, perturbations cor- 
responding to a tn as short as 10 4 £d may not suffice to yield 
a convergence of moments for two different chaotic ensem- 
bles with the same 7 and E. This is, e.g., illustrated in FIG- 
URE 12, which exhibits a z for the same two sets of initial 
conditions used to generate FIGURE 7, but now allowing 
for friction and white noise with tn — 10 to, 10 5 to, and 
10 4 to- The t — 200ir> dispersions for the two ensembles are 
appreciably different even for tn = 10 4 to, where the fric- 
tion and noise have caused the energy of a typical orbit to 
change by more than 4%. Both ensembles began with ener- 
gies Ex— 0.996. When evolved allowing for perturbations 
corresponding to tn — 10 4 io, the two final energy disper- 
sions were <te = 0.045 and erg = 0.046. 
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Figure 11. (a) The dispersion u z computed for the same initial 
conditions as FIG. 8, now allowing for coloured noise with = 
-E, ij = 0.4x 10 -5 , and t c = 1/a = 1.0 « 0At D . (b) The same 
for n = 0.4 X 10~ 5 , and t c = 3.0 1.2t D . (c) The same for 
v = 0.4 X 10~ 5 , and t c = 10.0 4t D . (d) The same for r) = 
0.4x 10~ 4 , and t c = 10.0 a4t D . (c) The same for n = 0.4 X 10~ 3 , 
and t c = 10.0si4t D . (f) The same for r) = 1.2 X 10~ 3 , and 
( c = 30.0 12*d. 



Figure 12. (a) The dispersion c z computed for one initially lo- 
calised ensemble of chaotic orbits in the lowest energy shell with 
7 = 1, allowing for friction and white noise corresponding to 
t/j = lQfttjj (b) a z for another chaotic ensemble with the same 
7, E and t R . (c) and (d) c z for the same two ensembles of ini- 
tial conditions, now integrated with tji = 10 5 t_o.(e) and (f) a z 
for the same two ensembles of initial conditions, integrated with 
t R = 10 4 t D . 



5 THE EFFECTS OF A SUPERMASSIVE 
BLACK HOLE 

The objective here is to determine how the basic picture 
described in the preceding Sections is altered if one allows 
for a modified potential of the form 



V(r) = V 7 (r) 



Me 



(x 2 +y 2 + z 2 + e 2 y/ 2 ' 



(10) 



where Vj denotes the unperturbed potential associated with 
the density distribution (1). Mbh is the black hole mass, 
expressed in units of the total galactic mass (since eq. (1) 
implies that J d 3 r p(r) = 1) and e = 1CP 3 is a soften- 
ing parameter. The choice of a relatively large value for e 
was motivated by computational considerations. Noisy sim- 
ulations are effected with a fixed time step integrator (cf. 
Honerkamp 1994), but integrations with a force which be- 
comes very large for r — > require very small time steps 
to maintain reasonable accuracy. Indeed, even for this com- 
paratively large value of e, the black hole simulations with 
Mbh = 1CP 3 and 7=1 required a time step an order 
of magnitude smaller than for comparable simulations with 
Mbh = 0. An investigation of the effects of varying e in the 
simpler potential (11), which exhibits qualitatively similar 
behaviour, indicates that the precise value of e is compar- 
atively unimportant, provided only that r min , the closest 
approach of an orbit to the black hole, is large compared 



with e. This condition suggests that it suffices to consider 
e< 1(T 2 . 

Four sorts of experiments were performed, namely: (1) 
computations of "libraries" of 1000 representative orbits of 
length t — 200tj3, to determine the relative numbers of reg- 
ular and chaotic orbits; (2) long time integrations for t = 
102400tn to obtain estimates of the true x fo r chaotic orbits, 
as well as to extract distributions of short time Lyapunov 
exponents; (3) 1600-orbit simulations of chaotic mixing, al- 
lowing for a black hole but no additional perturbations; and 
(4) 1600-orbit simulations of chaotic mixing that included 
both a black hole and low amplitude time-dependent irreg- 
ularities. Attention focused primarily on two choices of black 
hole mass, namely Mbh = 10~ 3 and Mbh = 10~ 2 . 



5.1 How chaotic are the orbits? 

Consistent with Merritt (1998), it was found that allowing 
for a nonzero Mbh makes the flow 'more chaotic' in the 
sense that, for chaotic orbits, the typical values of the Lya- 
punov exponents \ increase. Not surprisingly, this effect is 
most pronounced in the lower energy shells where the or- 
bits tend to pass the closest to the hole. For all four val- 
ues of 7, introducing a black hole with a mass as large as 
Mbh = 10~ 2 only increases the typical values of \ in tne 
ninth shell by less than a factor of two. However, especially 
for smaller values of 7, the black hole can have much larger 
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Figure 13. (a) Short time x(* = 200i D ) for 1000 orbits in the 
lowest energy shell with 7 = 0, arranged by increasing values of 
X- From bottom to top, the three curves represent Mbh = 0, 
M BH = 10~ 3 , and M BH = 10~ 2 . (b) The same for the lowest 
energy shell with 7 = 0.3. (c) The second lowest energy shell with 
7 = 0. (d) The second lowest energy shell with 7 = 0.3. (e) The 
ninth lowest energy shell with 7 = 0. (f) The ninth lowest energy 
shell with 7 = 0.3. 

effects on the lowest two energy shells. For 7 = 0, a black 
hole with Mbh = 1CF 2 increases the values of \ in the low- 
est energy shell by an order of magnitude. For 7 = 1, the 
increase is by a factor of three. For 7 = 2, a black hole with 
mass Mbh = 1CP 2 increases the typical \ by less than a fac- 
tor of two, and Mbh = 1CP 3 has a comparatively minimal 
effect. 

The presence of a supermassive black hole can also 
impact significantly the relative abundance of "strongly 
chaotic" orbit segments with x significantly larger than zero. 
For the two less cuspy models, 7 = and 7 = 0.3, even a 
black hole mass as small as Mbh = 1CP 3 suffices to trig- 
ger a significant increase in the relative number of strongly 
chaotic segments, particularly for the lower energy shells. 
For 7 = 1 and At = 200tc, such a black hole increases 
significantly the number of chaotic orbits in the two low- 
est energy shells, but its effect on the ninth energy shell is 
much less pronounced. For 7 = 2, even a black hole as large 
as Mbh = 1CP 2 seems to have no appreciable effect on the 
relative number of strongly chaotic orbits in any of the three 
energy shells. 

The net conclusion, in agreement with Merritt (1998, 
1999), is that, under appropriate circumstances, supermas- 
sive black holes with Mbh ~ 10 -3 — 1CF 2 can significantly 
impact both the relative number of chaotic orbits that exist 
and the typical values of their largest short time Lyapunov 
exponent. Merritt 's argument that the principal effect of the 
supermassive black hole is to enhance the strength of the 



central cusp would seem to explain why, for 7 = 2, a black 
hole of fixed mass has a much smaller effect than for smaller 
values of 7: In this case, a very pronounced density cusp is 
already there! 

But what about the question of 'stickiness'? Does the 
presence of a supermassive black hole tend to make distinc- 
tions between regular and chaotic orbits more evident, or is 
it still possible for a chaotic orbit to 'look regular' for very 
long times? Here the first obvious point is that, viewed over 
times ~ lOOto — 200to, the computed values of x f° r an 
ensemble containing both regular and chaotic orbits do not 
divide neatly into two disjoint classes with small and large 
values of \ and nothing in between. 

This is, e.g., manifest from FIGURES 13 and 14, which 
summarise information about the numbers of orbit seg- 
ments with different values of x f° r integrations of duration 
t = 200tD- These FIGURES encapsulate information about 
the short time Lyapunov exponents computed for 36 differ- 
ent ensembles of 1000 initial conditions - one each for the 
first, second, and ninth energy shells for the four different 
values of 7 with Mbh = 0, Mbh = 10~ 3 , and Mbh = 10~ 2 . 
For each ensemble, the orbits were ordered by the size of the 
computed value of x and the resulting x's were then plotted 
as a function of particle number from 1 to 1000. Inspection 
of these FIGURES thus allows one to infer the overall range 
of x's that was observed, the relative number of orbits in any 
given interval Ax, and how this number varies with changes 
in E, Mbh, and/or 7. In interpreting these FIGURES, it 
should be noted that they were generated as collections of 
tiny diamonds, not solid lines. The near-complete absence of 
any gaps thus manifests the fact that the orbits seemingly 
sample a near-continuous range of x' s , an d that there is no 
pronounced break between regular and chaotic orbits. The 
presumption is that all points not lying along a nearly hor- 
izontal line with very small x correspond to chaotic orbits, 
some relatively sticky and some otherwise. 

The absence of a pronounced gap in x between the reg- 
ular and chaotic orbits indicates that, as for the case with 
Mbh = 0, the distribution of short time x's must extend 
down to zero and, as such, one might perhaps suppose that 
JV[x], the distribution of short time Lyapunov exponents for 
chaotic orbit segments, is again at least bimodal. That this is 
often so is illustrated in FIGURES 15 and 16, which exhibit 
the analogues of several panels from FIGURES 1 and 2 for 
long time integrations with Mbh = 10~ 3 and Mbh = 10~ 2 . 

But how do things vary with the sampling time At? 
Does the presence of a supermassive black hole serve to ac- 
celerate phase space transport significantly, thus making the 
time scale to diffuse from one portion of phase space to an- 
other much shorter? The answer here seems to be: no! One 
again discovers that, over very long time scales, t ~ 20000t.o 
or longer, different segments of the same chaotic orbit can 
look extremely different, varying in appearance from wildly 
chaotic to nearly regular. This is manifested by the time- 
dependence of the dispersions a x (At) which, as for the case 
when Mbh = 0, typically decay much more slowly than as 
At^ 1 / 2 . Examples of the observed behaviour, generated as 
for FIGURES 3 and 4, are exhibited in FIGURE 17, which 
show a x for several different choices of 7, E, and Mbh- 
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Figure 14. (a) Short time x(* = 200*1?) for 1000 orbits in the 
lowest energy shell with 7 = 1. From bottom to top, the three 
curves represent M BH = 0, M BH = 10~ a , and M BH = 10~ 2 . (b) 
The same for the lowest energy shell with 7 = 2. (c) The second 
lowest energy shellwith7 = 1. (d) The second lowest energy shell 
with 7 = 2. (e) The ninth lowest energy shell with 7 = 1. (f) The 
ninth lowest energy shell with 7 = 2. 



5.2 Chaotic mixing 

The effects of a supermassive black hole on chaotic mixing 
are completely predictable. The presence of such a black 
hole fails to facilitate significantly the evolution of initially 
localised ensembles towards a true equilibrium in the sense 
that the phase space is still very sticky, so that chaotic or- 
bits can remain confined in restricted phase space regions 
for very long times. However, a supermassive black hole does 
tend to accelerate the approach towards a near-equilibrium. 
The correlation between the Lyapunov time th = 1/x an d 
the time scale tu on which a chaotic ensemble evolves per- 
sists for Mbh 7^ 0, but this means that larger x's correspond 
to a more rapid evolution. One example of this behaviour is 
shown in FIGURE 18, which exhibits a x for the same ini- 
tial ensembles evolved with 7 = 2 and variable Mbh = 0, 
Mbh = 1CT 3 , and Mbh = 1CT 2 . 

The effects of weak perturbations on chaotic mixing 
with a supermassive black hole are also predictable. As is 
the case for Mbh = 0, such perturbations can facilitate 
phase space transport, thus accelerating the approach to- 
wards equilibrium; and, even in the absence of any obvious 
'stickiness,' they can play an important role in 'fuzzing out' 
high frequency structures. Moreover, the minimum ampli- 
tudes required to have an appreciable effect are compara- 
ble to what is required when Mbh = 0. In particular, dis- 
creteness effects, modeled as friction and white noise, should 
be important on a time scale < lOOtc provided only that 
t R < 10 6 - 10 7 t D . 
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Figure 15. (a) N[x(At)] , the distribution of short time Lyapunov 
exponents, for chaotic orbits in the lowest energy shell with 7 = 
and Mbh = 10~ 3 for sampling time At = 100t o . (b) N[x{A)] for 
chaotic orbits in the lowest energy shell with 7 = 0, Mbh = 10~ 2 , 
and At = lOOt/j. (c) The same for chaotic orbits in the ninth 
energy shell with 7 = and Mbh = 10~ a (d) The same for 
chaotic orbits in the ninth energy shell with 7 = and Mbh = 
10~ 2 (c) - (h) The same as (a) - (d) for 7 = 0.3. 



6 DISCUSSION 

The 'stickiness' observed in these triaxial Dehnen potentials 
is so striking that one may wonder how generic it is. It is thus 
reassuring that the same qualitative behaviour can arise for 
the simple model comprised of an anisotropic oscillator and 
a spherical Plummer potential. For example, the qualitative 
form of iV[x(A£)], the distribution of short time Lyapunov 
exponents for the triaxial Dehnen potentials, can be repro- 
duced by orbits evolved in a potential of the form 

Mbh 

■ Jy y - ■ 



V(r) 



22, 22. 22 

Ur.x +u y y +oj z z 



Vr 2 + e 2 



(11) 



with to 2 = 1.0, uo 2 y = 1.25, to' 2 = 0.75, and e = 10~ 3 , for 
appropriate choices of energy E and black hole mass Mbh- 
For Mbh = and Mbh — * 00, motion in this poten- 
tial is completely integrable. However, V(r) is nonintegrable 
for intermediate values and, for a reasonable range of black 
hole masses, admits a coexistence of large numbers of both 
regular and chaotic orbits. In this range, one finds typically 
that increasing Mbh makes the system more chaotic in the 
sense that the values of \ tend to increase, but that, when 
expressed in units of l/t D , the changes in x are n °t all that 
large. One discovers further that, generically, the chaotic or- 
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Figure 16. (a) N[x(At)], the distribution of short time Lyapunov 
exponents, for chaotic orbits in the lowest energy shell with 7 = 1 
and M BH = 1(T 3 for sampling time At = 100t D . (b) N[x(A)] for 
chaotic orbits in the lowest energy shell with 7 = 1, Mbh = 10~ 2 , 
and At = lOOt/j. (c) The same for chaotic orbits in the ninth 
energy shell with 7 = 1 and Mbh = 10~ 3 (d) The same for 
chaotic orbits in the ninth energy shell with 7 = 1 and Mbh = 
1CT 2 (c) - (h) The same as (a) - (d) for 7 = 2. 




Figure 17. (a) The dispersion a x for chaotic orbits in the lowest 
energy shell with 7 = 1 and Mbh = 10~ 3 , computed as a func- 
tion of sampling time At. (b) a x for chaotic orbits in the lowest 
energy shell with 7 = 1 and Mbh = 10~ 2 . (c) The same for the 
ninth energy shell and Mbh = 10~ 3 . (d) The same for the ninth 
energy shell and Mbh = 10~ 2 . (c) - (h) The same as (a) - (d) 
for 7 = 2. 



bits tend to be very sticky, so that an integration for a time 
t — lOOOOin or longer does not suffice to yield convergence 
towards the asymptotic x, defined in a t —> oo limit. In both 
these ways, this potential is very Dehnenesque. 

As an indication of the degree to which this poten- 
tial yields similar distributions of short time Lyapunov ex- 
ponents N[x{At)], consider FIGURE 19, which exhibits 
N[x(At)] for ensembles of chaotic orbits. These ensembles 
were generated from initial conditions that uniformly sam- 
pled the y — z plane, setting x — v v = v z = and, for given 
energy E, solving for v x (y,z,E) > 0. Each orbit was inte- 
grated for a time t = 4096 and the computed values of \ 
at that time were used to identify which initial conditions 
appeared chaotic. The orbits identified as chaotic were then 
divided into 16 segments of length At = 256 which, for the 
range of energies and black hole masses exhibited in these 
FIGURES, corresponded to a period of order lOOto, and 
short time Lyapunov exponents were identified for each seg- 
ment. The resulting collection of x's was then binned to gen- 
erate N[x] for At = 256. The first panel, for Mbh = 10" l s 
and E — 0.0044, closely resembles the bimodal distributions 
iV[x(A£)] observed for the lower energy shells in the Dehnen 
potentials. The second panel, generated for Mbh = 10~ 2 



and E = 0.65, resembles more closely the distributions for 
the higher energy shells, especially in the presence of a black 
hole. 

The results described in this paper have poten- 
tially significant implications for modeling galaxies using 
Schwarzschild's (1979) method or any comparable tech- 
nique. One obvious point is that, in the absence of any time- 
dependent perturbations, phase space transport can be ex- 
tremely inefficient, so that computing an ensemble of orbits 
even for lOOOtn or longer may not suffice to obtain a rea- 
sonable approximation to an invariant distribution. Great 
care needs to be taken in endeavoring to create truly time- 
independent building blocks. At least for cusps that are 
not too steep, inserting a supermassive black hole typically 
makes orbits more chaotic in the sense that the Lyapunov 
exponents become larger and, in addition, tends to increase 
the relative abundance of chaotic orbits. However, the pres- 
ence of the black hole does not seem to significantly accel- 
erate the rate of phase space transport. 

The other obvious point is that even comparatively 
weak time- dependent perturbations tend to significantly ac- 
celerate phase space diffusion, allowing orbits to probe dif- 
ferent phase space regions on times short compared with the 
time scale for phase space transport in the absence of any 



© 1999 RAS, MNRAS 000, 000-000 



16 C. Siopis and H. E. Kandrup 



0.030 
0.020 

0.010 
0.000 




(a) 







25 

t/t n 



50 



0.030 
0.020 

0.010 
0.000 




(c) 







25 
t/t D 



50 



0.030 
0.020 

0.010 
0.000 




(e) 







25 
t/t„ 



50 




Figure 18. (a) The dispersion a x computed for an unperturbed 
ensemble of 1600 orbits in the second lowest energy shell for 7 = 2, 
evolved in the absence of a supcrmassivc black hole, (b) The same 
for an ensemble in the ninth energy shell, (c) The same as (a), 
but now allowing for a black hole with mass M BH = 10- 3 . (d) 
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Figure 19. (a) Short time N[x(t = 256)] for a collection of 
chaotic orbits evolved in the toy potential (11) with Mbh = 
10- 18 and E = 0.044. (b) The same for M BH = 10~ 2 and 
E = 0.065. 



perturbations. This basic fact, first recognised in the con- 
text of simple maps (cf. Lieberman & Lichtenberg 1972) and 
subsequently explored for simple two- and three-degree-of- 
freedom systems (cf. Habib, Kandrup, & Mahon 1997), per- 
sists for seemingly realistic cuspy, triaxial potentials. This 
suggests immediately that irregularities associated with the 
real world could serve to destabilise pseudo-equilibria con- 
structed using near- invariant, rather than truly invariant, 
building blocks, even if, in the absence of perturbations, such 
equilibria could persist for times 3> tn- 

But what are the effects of perturbations like friction 
and noise? Do they simply accelerate the evolution of an 
unperturbed system or do they alter the final state? Noise in- 
duces significant changes in the energy and any other global 
isolating integrals on a relaxation time tR = 1/r) and, as 
such, it is clear that noise must change the form of the final 
state for later times. However, there is good reason to antic- 



ipate that, for times t <§C tR, noise simply speeds things up 
without significantly altering the ultimate outcome. In other 
words, an ensemble evolved for a time tn in the presence of 
weak noise with tR >• tn might be expected to achieve a fi- 
nal state closely resembling the state which an unperturbed 
ensemble would have achieved only after a time lOOiff, or 
even longer. The basic idea is that the perturbations drive 
phase space transport by helping stars to find phase space 
holes, thus accelerating the approach towards an equilib- 
rium, but that, at least to the extent that the perturbations 
do not significantly alter the values of the energy and other 
isolating integrals, they will not impact the form of the final 
equilibrium. 

An analogy with kinetic theory were perhaps appro- 
priate. When analysing interactions of particles in a dilute 
gas, one speaks of microreversibility, which implies that the 
transition rates W(vi, V2 — + V3, V4) and W(\ 3 , V4 — > vi, V2) 
must be equal. This is a statement about the phase space 
kinematics, and is without reference to the dynamics. Sta- 
tistical equilibrium requires detailed balance, characterised 
by populations /(v<) satisfying /(vi)/(v2)W(vi, V2 — > 

V 3 ,V 4 ) = /(V 3 )/(V 4 )W(V3,V4 -> Vi,V 2 ). 

In the context of phase space diffusion, the analogue 
of microreversibility is the notion that the transition rates 
W(A —> B) and W(B — > A) across a leaky barrier sepa- 
rating regions A and B must be equal. Equilibrium involves 
a detailed balance with equal (and constant) number den- 
sities on both sides of the boundary, so that the product 
of number density and transition rate is the same. That 
W(A — + B) = W(B — > A) is incorporated explicitly in every 
theory of phase space transport, including, e.g., the "turn- 
stile model" for diffusion through cantori, which is the most 
successful to date (MacKay, Meiss, & Percival 1984). The 
expectation, then, is that low amplitude perturbations will 
help stars breach these leaky barriers, but that they will 
not alter the balance associated with an equilibrium. One 
could, perhaps, envision carefully tailored "designer noise" 
so constructed as to favour motion in one direction over the 
other, but this would seem contrived. The noise used in this 
paper, and most other investigations of accelerated phase 
space transport, was chosen to act equally on stars on both 
sides of the boundary; and, as such, should not serve to alter 
the final balance. 

There is also concrete numerical evidence (cf. Habib, 
Kandrup, & Mahon 1996, 1997) suggesting that weak noise 
simply accelerates phase space transport rather than leading 
to a new equilibrium. This derives from comparing the same 
ensemble of initial conditions evolved in three different ways, 
viz: (i) for comparatively short times, t < lOOto, in the ab- 
sence of any perturbations, (ii) for longer times, t > lOOOtrj, 
in the absence of perturbations, and (iii) for short times, 
again t < WOto, in the presence of weak noise for which 
tR 3> lOOtrj. For some choices of potential and/or ensemble, 
one finds that the short and long time unperturbed simu- 
lations result in very similar end states, as probed, e.g., by 
lower order moments and coarse-grained distribution func- 
tions. In this case, weak noise has only a minimal effect. 
In other cases, however, the short and long times unper- 
turbed simulations yield substantially different end states. 
In these cases, one finds invariably that the noisy simula- 
tions yield a final state that more closely resembles the late 
time unperturbed state than did the shorter time unper- 
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turbcd evolution. Moreover, for noise of sufficient strength, 
one discovers oftentimes that the noisy end state and the 
long time unperturbed endstate are almost identical. 

This motivates a specific proposal, namely that, when 
using Schwarzschild's method to construct models that con- 
tain significant numbers of chaotic orbits, one should gen- 
erate "noisy libraries," which are more likely to constitute 
reasonable approximations to invariant distributions than 
are libraries generated without any perturbations. There is 
no guarantee that the noise will yield truly time-independent 
building blocks. However, it seems reasonable to anticipate 
that allowing for weak noise will yield libraries that are, at 
least, better approximations to invariant distributions. 

In contemplating the implications of the computations 
described in this paper for real galaxies, it is important 
to recall that the notion of an equilibrium is an idealisa- 
tion which is never achieved in the real world. Assuming 
that the evolution is not strongly dissipative, so that Li- 
ouville's Theorem holds at least approximately, there can 
at best be a coarse-grained approach towards an equilib- 
rium (cf. Lynden-Bell 1967). One might anticipate an ap- 
proach towards equilibrium in the sense that various mo- 
ments and/or coarse-grained distributions become system- 
atically more nearly time-independent, but there can be no 
true pointwise approach towards equilibrium. Moreover, re- 
alistic systems are continually subjected to various time- 
dependent perturbations reflecting both internal and exter- 
nal irregularities, so that, even allowing for dissipation, one 
could at best view a galaxy as being "near" equilibrium. The 
basic question is whether it be reasonable to assume that a 
real galaxy evolves towards a complex time-dependent state 
which can be approximated realistically as an equilibrium. 

Evolving towards a true triaxial equilibrium might seem 
relatively difficult! It would seem unreasonable to model re- 
alistic triaxial systems as equilibria that depend only on 
one global integral and nothing else. A one-integral f{E) 
depending only on energy E must be spherical (cf. Perez & 
Aly 1996). It is possible to construct rotating triaxial equi- 
libria f(Ej) which involve only the Jacobi integral Ej, but 
such equilibria apparently cannot be sufficiently centrally 
condensed to mimic real physical systems (cf. Vandervoort 
1980, Ipser & Managan 1981). Potentials that admit two or 
three global integrals are a set of measure zero in the set 
of three-dimensional potentials, and there is a precise sense 
in which near-integrable systems are substantially rarer for 
three- and higher-degree-of-freedom systems than for sys- 
tems with two degrees of freedom (Poincare 1892). 

There remains the possibility of trying to construct 
equilibria involving non-standard local integrals in addition 
to the conventional global isolating integrals, perhaps ex- 
cluding the chaotic orbits altogether. This possibility seems 
implicit in virtually all work on cuspy triaxial systems hith- 
erto (cf. Zhao 1996, Siopis 1998) and an attempt has been 
made recently to do this in a more systematic fashion 
(Hafher et al 1999). Most, if not all, of the potentials that 
have been considered, including the triaxial Dehnen poten- 
tials, have chaotic orbits with two positive Lyapunov expo- 
nents and, consequently, can admit only one global isolating 
integral. However, many dynamicists seem to have the in- 
tuitive expectation that such equilibria would be hard to 



realise, since they entail a "fine-tuning" on a comparatively 
microscopic level. As such, one might therefore ask: is it rea- 
sonable to expect that the true time-dependent distribution 
function associated with a real galaxy actually manages to 
approach one of these finely tuned distributions involving 
local integrals? 

Irrespective of the answer to this question, it would 
seem substantially easier to approach a near-equilibrium 
than a true equilibrium. One would anticipate intuitively 
that there are many more near-equilibria than true equilib- 
ria, and these could be supported approximately with two 
types of building blocks, namely (i) building blocks that are 
only approximately time-independent and/or (ii) building 
blocks which, albeit exactly time-independent, only repro- 
duce the potential approximately. In this regard, it should 
perhaps be stressed that there is no guarantee that, in the 
absence of dissipation, a real galaxy will evolve towards a 
time-independent equilibrium. One could, at least in princi- 
ple, envision an evolution towards a final state that entails 
undamped, finite amplitude oscillations about some time- 
independent equilibrium (cf. Louis & Gerhart 1988, Sridhar 
1989). 

These observations are related to two significant limi- 
tations implicit in Schwarzschild's method or any obvious 
variant thereof: (1) There is no guarantee that any solution 
constructed using Schwarzschild's method corresponds to a 
true equilibrium. It is possible that, strictly speaking, the 
assumed potential admits no true equilibria, and that the 
numerical algorithm has only constructed an approximate 
equilibrium (or worse). (2) The fact that, for some assumed 
potential, one's orbit library could not be used to generate 
a Schwarzschild equilibrium does not imply that an equilib- 
rium does not exist; and even if there is no equilibrium this 
does not preclude the possibility that there exist "nearby" 
potentials that correspond to a system which (at least on 
the average) only evolves comparatively slowly. 

The enormous stickiness manifested by chaotic orbits 
in (at least some) cuspy triaxial potentials can facilitate the 
construction of approximate equilibria. However, there is ev- 
ery reason to think that low amplitude perturbations asso- 
ciated with both internal and external irregularities could 
trigger nontrivial evolutionary effects on a time scale <C tn- 

From this slightly different chain of reasoning, one is 
led to a conclusion similar to that of Merritt (1996), namely 
that galaxies could exhibit a two-stage evolution, a rapid ap- 
proach towards a roughly, but not truly, time-independent 
equilibrium on a time scale ~ tn which is succeeded by a 
substantially slower evolution as the system "searches" for 
a true equilibrium. Using this line of reasoning exclusively, 
there is no obvious guarantee that this slower evolution will 
lead to configurations more nearly axisymmetric - one might 
instead see an evolution through a sequence of compara- 
bly triaxial equilibria. However, the apparent fact that very 
cuspy galaxies tend to be more nearly axisymmetric, at least 
in the center, could reflect in part the fact that, in these high 
density central regions, dissipation was sufficiently strong to 
allow the system to "develop a global conserved quantity." 
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